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Abstract 

The modern imaging techniques of Positron Emission Tomography and of Single 
Photon Emission Computed Tomography are not only two of the most important 
tools for studying the functional characteristics of the brain, but they now also play 
a vital role in several areas of clinical medicine, including neurology, oncology and 
cardiology. The basic mathematical problems associated with these techniques are 
the construction of the inverse of the Radon transform and of the inverse of the so 
called attenuated Radon transform respectively. We first show that, by employing 
mathematical techniques developed in the theory of nonlinear integrable equations, 
it is possible to obtain analytic formulas for these two inverse transforms. We then 
present algorithms for the numerical implementation of these analytic formulas, 
based on approximating the given data in terms of cubic splines. Several numerical 
tests are presented which suggest that our algorithms are capable of producing 
accurate reconstruction for realistic phantoms such as the well known Shepp-Logan 
phantom. 

1 Introduction 

Positron emission tomography (PET) and single photon emission computed tomogra- 
phy (SPECT) are two modern imaging techniques with a wide range of medical applica- 
tions. Although these techniques were originally developed for the study of the functional 
characteristics of the brain, they are now used in many diverse areas of clinical medicine. 
For example a recent editorial in the New England Journal of Medicine ^ emphasized 
the importance of PET in oncologic imaging. Other medical applications of PET and 
SPECT are presented in |2]-|22j. 

The first step in PET is to inject the patient with a dose of a suitable radiopharmaceuti- 
cal. For example in brain imaging a typical such radiopharmaceutical is flurodeoxy glucose 
(FDG), which is a normal molecule of glucose attached artificially to an atom of radioac- 
tive fluorine. The cells in the brain which are more active have a higher metabolism, need 
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more energy, thus will absorb more FDG. The fluorine atom in the FDG molecule suffers a 
radioactive decay, emitting a positron. When a positron collides with an electron it liber- 
ates energy in the form of two beams of gamma rays travelling in opposite direction, which 
are picked by the PET scanner. SPECT is similar to PET but the radiopharmaceuticals 
decay to emit a single photon. 

In both PET and SPECT the radiating sources are inside the body, and the aim is 
to determine the distribution g{xi,X2) of the relevant radiopharmaceutical from measure- 
ments made outside the body of the emitted radiation. If f{xi, X2) is the x-ray attenuation 
coefficient of the body, then it is straightforward to show [SB] that the intensity / outside 
the body measured by a detector which picks up only radiation along the straight line L 
is given by 

1 = J e~-^M-)^%dr (1.1) 

where r is a parameter along L, and L{x) denotes the section of L between the point 
{xi,X2) and the detector. The attenuation coefficient f{xi,X2) is precisely the function 
measured by the usual computed tomography. Thus the basic mathematical problem in 
SPECT is to determine the function g{xi,X2) from the knowledge of the "transmission" 
function f{xi,X2) (determined via computed tomography) and the "emission" function I 
(known from the measurements). 

In PET the situation is simpler. Indeed, since the sources eject particles pairwise in 
opposite directions and the radiation in opposite directions is measured simultaneously, 
equation ()1.1|) is replaced by 

1 = J e"^^+w^'^'"-^^-W^%dr, (1.2) 

where L+, L_ are the two half-lines of L with endpoint x. Since L+ + L_ = L, equation 
()1.2|1 becomes 

We recall that the line integral of the function /(xi, X2) along L is precisely what is known 
from the measurements in the usual computed tomography. Thus since both / and the 
integral of f{xi,X2) are known (from the measurements of SPECT and of computed to- 
mography respectively), the basic mathematical problem of PET is to determine g{xi, X2) 
from the knowledge of its line integrals. This mathematical problem is identical with the 
basic mathematical problem of computed tomography. 

Notation 

(i) A point of a line L making an angle 6 with the Xi-axis is specified by the three real 
numbers {t,p,6), where r is a parameter along L, —00 < r < 00, p is the distance from 
the origin to the line, —00 < p < 00, and < 6 < 2tt. 

(ii) The above parameterization implies that, for a fixed 6, the Cartesian coordinates 
(xi,X2) can be expressed in terms of the local coordinates {t, p) by the equations (see 
Section 

xi = T cos6 — psinO, X2 = t sin6 + pcosO. (1-3) 
A function f{xi,X2) rewritten in local coordinates will be denoted by F{t,p,6), 
F{t, p, 6) = f{r cos 6* — psin 6*, r sin 9 + pcos 9). 
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Figure 1: Local coordinates for the mathematical formulation of PET and SPECT. 



Thus F{t,p,6) and G{t,p,6) will denote the x-ray attenuation coefficient f{xi,X2) and 
the distribution of the radiopharmaceutical g{xi,X2), rewritten in local coordinates, 
(iii) The line integral of a function / is called its Radon transform and will be denoted by 
/. In order to compute /, we first write / in local coordinates and then integrate with 
respect to r, 

/■oo 

fip,0)= / FiT,p,e)dT. (1.4) 



The line integral of the function g with respect to the weight / appearing in equation (jl.lj) 
is called the attenuated Radon transform of g (with the attenuation specified by /) and 
will be denoted hj gj. In order to compute gf,^e write both g and / in local coordinates 
and then evaluate the following integral 

/oo 
^-rF(.,p,.)d.^(^^^^^)^^_ (1.5) 
-oo 

Mathematical Methods 

The basic mathematical problem of both computed tomography and PET is to recon- 
struct a function / from the knowledge of its Radon transform /, i.e. to solve equation 
()1.4|) for f{xi,X2) in terms of f{p,6). The relevant formula is called the inverse Radon 
transform and is given by 

/(^i,^2) = -^(9., -i5.J Te'^ (f r^^^4^^—^] (1-6) 

where — oo < Xj < oo, j = 1, 2 and ^ denotes principal value integral. 

A novel approach for deriving equation (jl.fij) was introduced in j21], and is based on 
the analysis of the equation 

\ + \] dx^ + ^ ( X-\] d^A K^i^^-2,>^) = f{xi,X2), (1.7) 
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where A is a complex parameter different than zero. The application of this approach to 
a slight generalization of equation ()1.7p can be used to reconstruct a function g from the 
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knowledge of its attenuated Radon transform gj, i.e. this approach can be used to solve 
equation ()1.5|) for g{xi, X2) in terms of gf{p, 0) and /(xi, X2). The relevant formula, called 
the inverse attenuated Radon transform, was obtained by R. Novikov jSH] by analysing, 
instead of equation p.7|l . the equation 

\ + \]d:,^ + ^{\-\]d^^ + f{xi,X2)]p.{xi,X2,\) = g{xi,X2). (1.8) 



2 V a; ^ 2i V 

Organization of the Paper 

In Section|21we first review the analysis of equation ()1.7|1 . and then show that if one uses 
the basic result obtained in this analysis, it is possible to construct immediately the inverse 
attenuated Radon transform. In Section El we present a new numerical reconstruction 
algorithm for both PET and SPECT. This algorithm is based on approximating the given 
data in terms of cubic splines. We recall that both the exact inverse Radon transform 
as well as the exact inverse attenuated Radon transform involve the Hilbert transform of 
the data functions. For example, the inverse Radon transform involves the function 

h{p,e)= f l^dp'. (1.9) 

/_oo P'-P 

Existing numerical approaches use the convolution property of the Fourier transform to 
compute the Hilbert transform and employ appropriate filters to eliminate high frequen- 
cies. It appears that our approach has the advantage of simplifying considerably the 
mathematical formulas associated with these techniques. Furthermore, accurate recon- 
struction is achieved, for noiseless data, with the additional use of an averaging or of 
a median filter. Several numerical tests are presented in Section HI One of these tests 
involves the Shepp-Logan phantom |^, see Figure Efc). 

Numerical algorithms based on the filtered back projection are discussed in pTj-pU], 
while algorithms based on iterative techniques can be found in 



2 Mathematical Methods 

We first review the basic result of ||24j. It will be shown later that using this result 
it is possible to derive both the inverse Radon as well as the inverse attenuated Radon 
transforms in a straightforward manner. 

Proposition 2.1. Define the complex variable z by 

where xi, X2 are the real Cartesian coordinates —00 < Xj < 00, j = 1,2, and A is 
a complex variable, A 7^ 0. Assume that the function f{xi,X2) has sufficient decay as 
+ \x2\ —>■ 00. Let fi{xi,X2,X) satisfy the equation 

^(^-|Ar)^^^^^||^ = /(x„X2), |A|^1, (x„X2)GM^ (2.2) 
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as well as the boundary condition /i = 0(1/2;) as |a;i| + |x2| — > oo. Let A"*" and A~ 
denote the limits of A as it approaches the unit circle from inside and outside the unit 
disc respectively, i.e. 

\^ = \im{l T e)e'\ £ > 0, < ^ < 27r. 

e—*0 




Figure 2: The unit circle. 

Then 

/oo 
F{T',p,9)dT', (2.3) 

where / denotes the Radon transform of /, F denotes / in the local coordinates (see the 
Notation in Section H}, P"^ denote the usual projection operators in the variable p, i.e. 

^ j-o2ri/_„p'-(p±i£) 2 2mf_^p'-p' ^ ' 

and ^ denotes the principal value integral. 

Proof. Before deriving this result, we first note that equation ()2.1|) is a direct consequence 
of equation p.7|) . Indeed, equation ()1.7p motivates the introduction of the variable z 
defined by equation ()2.1|) . Taking the complex conjugate of equation ()2.1|) we find 



Equations 1)2.11) and ()2.5|1 define a change of variables from (xi,X2) to {z,z). Using this 
change of variables to compute d^^ and in terms of dz and d^, equation p.7|) becomes 



We now derive equation ()2.3)1 . The derivation is based on the following two steps, which 
have been used extensively in the field of nonlinear integrable PDEs, see for example [34 . 

(i) In the first step (sometimes called the direct problem), we consider equation ()2.2j) 
as an equation which defines /i in terms of /, and we construct an integral representation 
of fj, in terms of /, for all complex values of A. This representation is 

pix^,x,,X) = ^sgn(^^ - |Ap^ J I' M^dx[dx',, |A| ^ 1. (2.6) 

R2 

Indeed, suppose that the function p{zr, zj) satisfies the equation 

dp{zR, zi] 



dz 



g{zR,zi), z = zr + izi, -oo<zr<(X), -oo < zj < oo, 
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as well as the boundary condition /i = 0(1/^;) as ^ — >■ oo. Then Pompieu's formula (see 
for example j^) implies 

W|9(£^ (2.7) 

nJJ z'-z 



In our case 



^^■^ ^zr^zi = f TTT^ - ) dxida;2, 



^-|A|2' " ' 2iVlA 



thus equation ()2.7p becomes ()2.6|) . 

(ii) In the second step (sometimes called the inverse problem), we analyze the ana- 
lyticity properties of with respect to A, and we find an alternative representation for 
H. This representation involves certain integrals of / called spectral functions. For our 
problem, this representation is equation ()2.3p . Indeed, since /i is an analytic function of 
A for |A| 7^ 1 and since fi = 0(1/A) as A ^ oo, we can reconstruct the function fi if we 
know its "jump" across the unit circle: 

^{xi,X2,X) = — ^.g, _ ^ d^, (2.8) 

where 

J{xi,X2,9) = /i(a;i,X2, A+) - fi{xi,X2, X~). 
Thus we need to compute the limits of /x as A tends to A^. As e — > 0, 

Substituting this expression in the definition of z (equation 1)2.11) ) and simplifying, we find 



z' — z {x'l — Xi) sin 6* — (x'a — X2) cos9 + ie{{x[ — Xi) cos 6* + (^2 — X2) sin 6*). (2.9) 

The right-hand side of this equation can be rewritten in terms of the local coordinates p, 
p', r, r': Let k and k"*- denote two unit vectors along the line L and perpendicular to this 
line, respectively. Then 

X = rk + pk"*", 

or 

(xi, X2) = r(cos 9, sin 9) + p(— sin 9, cos 6*). 
Hence xi and X2 are given by equations ()1.3)1 . Inverting these equations we find 

r = X2 sin6' + xi cos6', p = X2 cos6' — xi sin6'. (2.10) 

Thus equation ()2.9p becomes 

z' ~ z —p + p + ieir' — r). 

Substituting this expression in equation ()2.6p and using the fact that the relevant sign 
equals 1, we find 

Mxi,X2,A+)^--L // ^_^o^ ^^o_ (2.11) 

2-Ki J J p' - p - ie[T' - t) 
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Using the change of variables (xi,X2) ^ (t, p) defined by equations ()1.3p and ()2.10p . and 
noting that the relevant Jacobian is 1, i.e. 

f{x[,x',)dx[dx', = F{T',p',9)dT'dp', 

we find that the right-hand side of equation ()2.11|1 equals 

1 rr fdrv 



2m J J p' - {p + ie(r' - r)) ' 
In order to simplify this expression we split the integral over dr' in the form 

/»r poo 

dr' = / dr'+ / dr', 



and note that in the first integral t' — t < 0, while in the second integral t' — t > 0. 
Thus, using the second set of equations ()2.4|) the expression in ()2.12|) becomes 

Finally, adding and subtracting the integral | we find 

Mxi,X2,A+) = r ff FiT',p',9)4^)dT' 

27ri y_oo V/-00 P - py 



-1 poo poo 

+ - j F(r',p,^)dr'- j F(r',p,^)dr'. 



The first two terms in the right-hand side of this equation equal —P /, hence we find 
(12.31) "*^. The derivation of equation ()2.3p ~ is similar. QED 

Using equation ()2.3|1 it is now straightforward to derive both the inverse Radon and 
the inverse attenuated Radon transforms. In this respect we note that the result of 
Proposition 12.11 can be rewritten in the form 

Hm {^J^ (^-^^^y^=TP^f{p.e)~ F{^^ (2.13) 



where 



KA) = ^(^-|An. (2.14) 



2i V A 



The Inverse Radon Transform 
Equations ()2.3|) yield 



j(..,..,e).-^f—Mf^^. (2.15, 

^1 J -00 P ~ COS U — Xi sm 6) 



Equation ()2.8j) implies 



1 .„\ 1 „ / 1 



p(a;i, X2, A) = ( -— / J(xi, X2, 9)e''d9 ) ^ + O ( - 
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Substituting this expression in equation ()1.7p we find 



(2.16) 



Replacing in this equation J by the right-hand side of equation ()2.15|) we find equation 
(Hil). 



The Attenuated Radon Transform 

Equation p.8j) can be rewritten in the form 

dfi f g 

OZ V V 

where v is defined by equation ()2.14j) . Hence 



d 

-,/.exp 



V 



9 

— exp 

V 



or 



/xexp 



V 



V 



exp 



Replacing in this equation 9^ ^ (^) by the right-hand side of equation ()2.13j) we find 

V J 

For the computation of the right-hand side of this equation we use again equation ()2.13|) . 
where / is replaced by g times the two exponentials appearing in the above relation. 
Hence 

/.(Xi, X2, A±)e^^"^"(^'^)e- ^ ^r'^P^r' = 

POO 

^pT^TP^f(pfi)^g^^p^e)- / G(r'p,^)e^^^^"(''''')e-/^°^^(^''''^)"Mr'. (2.17) 



Note that the term exp[=pP'^/] is independent of r', thus this term comes out of the 
integral J°°, and furthermore the same term appears in the left-hand side of equation 
()2.17|l . Hence when computing the jump /u(xi, X2, A"*") — fi{xi,X2, A~), the second term in 
the right-hand side of equation ()2.17|1 cancels and we find that the relevant jump in now 
given by 

J{x,,X2,e) = -e^rnr',P,ew (^^p-np,e)p-^-p-fip,e)^^-p+f(P,e)p+^P+f(P,o)^ g^(^p^e) 

(2.18) 

where r and p are expressed in terms of Xi and X2 by equations ()2.1U|) . 

Equation ()2.8p is still valid, furthermore equation ()2.1(ij) is valid if / is replaced by g. 
Hence replacing in equation ()2.1(j|l f hj g we find 

r.27r 

J(xi,a;2,0)e'^d^, (2.19) 
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g[xuX2} 



where J is defined by equation ()2.18|1 . This formula is equivalent to Novikov's formula. 

In summary, letgf{p,9) be defined by equation / ti.,5|) . let F{t, p,9) denote the function 
f{xi,X2) written in local coordinates (see the Notation) and let f{p,6) denote the Radon 
transform of f{xi,X2) (see equation (|i.^| )j. Then g{xi,X2) is given by equation ^2.1fJ\) 
where the function J is explicitly given in terms of cjf and f by equation 
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3 Reconstruction Algorithm 

3.1 PET Algorithm 

Taking the real part of equation it follows that f{xi,X2) is given by 



f{Xl,X2) 



1 



(3.1) 



where h{p, 6) is defined by equation ()1.9|1 . 

We assume that /(xi, X2) has compact support, namely /(xi, X2) = 0, for xi'+X2^ > 1. 
For the numerical calculation of the integral in ()3.H) we use the formula 



27r 



27r 



N-l 



i=0 



27rz 



(3.2) 



Since g is analytic and periodic, this equispaced quadrature converges at spectral speed 
[HE] • In other words, ()3.2j) represents the optimal quadrature formula for the above integral 
and its implementation is likely to result in high precision even for relatively small values 
of A^. For the numerical calculation of /ip(p, 6) we suppose that /(p, 6) is given, for every 
6*, at n equally spaced points pi G [—1, 1], i.e. we suppose that fi = f{pi,0) are known. 
Moreover, in each interval [pi, pi^i] we approximate f{p,9) using the relation 



(3.3) 



where 



A = ^^±^^, B, = i-A, c, = ^(A,3-A,)(p,+l-p,)^ A = ^(5«'-5i)(p,,+i-A)', 

Pi+i - Pi 6 6 

and fl' denotes the second derivative of /(p, 9) with respect to p, at p = p^. In other words, 
we approximate f{p,9) by a cubic spline (in p) with equally-spaced nodes. Integrating 
the spline, we derive a well-known quadrature formula which, in our setting, reads 

n-l 



hp,9) = j2 / 



^^^^ S^{P\9)_^^, 



Following straightforward calculations we obtain 

n-l 



K{p, 9) 



i+l 



-M- 3p.+i + 2p)/;' - ^(3p, - p,+i - 2p)fl^ 

Pi- P Pi+i - p 4 4 



+ 



+ 



i=l 

fi ~ fi+1 _ 1 
Pi - Pi+i 6 
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Pi - Pi+i 



Pi - Pi+i 



Pi - Pi+1 



3(Pi+i - py 



n 



f" 



pi - pi+i 
pi+i - p 



+1 



In 



P^- P 



(3.4) 



In order to calculate numerically f{xi,X2) from the data f{p,9) we first compute the 
second derivatives /f. For this purpose we use the subroutine spline from Numerical 



9 



Recipes jHZI, setting f" = = (i.e. we use the natural cubic spline interpolation). 
Then, for any xi and X2, we calculate (for any 9) p using ()2.10b ) and hp{p, 6) using ()3.4|1 . 
Finally we calculate /(xi,X2) using (jH.lj) . 
We note that ()H.4|1 contains the term 



In 



Pi+i - P 



pi- P 



However, since for the reconstruction the number of the points for Xi and X2 can be 
different than the number of the p points, in general p 7^ p^+i and p ^ Pi- 



3.2 SPECT Algorithm 

We denote the first exponential term of the right-hand side of ()2.18|) by /(r, p, 9), i.e. 



P, = exp 



F{T',p,e)dT' 



(3.5) 



Note that, since we have assumed compact support, the integration domain is finite, i.e. 
[r, v^l and F(r, p, ^) = for |p| > 1, or for |r| > ^1 - p2. 

The definitions ()2.4|) become 

P±/(p,0) = ±^/(p,^^)--^%,0). 



Moreover 

exp P^f{p,6) =exp 
exp —P^f{p,9) = exp 
We introduce the following notation: 



±^fip,d) 



T^fiP, 



2tx 



hip,e) . . hip,e) 

cos 1 sm 

271 2tt 

h{p,9) , . . h{p,i 



cos ■ 



27T 



+ 1 sm 



27r 



(3.6) 



f-P^(p, 9) = e'^f^P^'^ cos /^^^(p, 0) = ei^'(^'^) sin 

/-'^(p, 9) = e-y^P''^ cos r™^(p, ^) = e-^^"(^'^) sin (3.7) 

27r 27r 

r(p, 9) = r%p, 9)gfip, 9), rip, 9) = r%p, 9)gjip, 9). (3.8) 

Using this notation and setting R{t,p,9) = — J{r,p,9), after some calculations, equation 
()2.18|1 becomes 

i?(r, p, 9) = I{t, p, 9) - if^^){p-r + iP-f) + + ir-^)(P+r - ^PT)) ■ 

(3.9) 

We now set 

/-oo P - P J-oc P -P 
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thus equation ()3.9p becomes 



R{t,p, 



i/(r, p, 9) /""^ -h' + 2/M + f'""' -h' - 2/ 



IT 



We denote the right-hand side of this equation by —ir{T,p,6). Taking the real part of 
g{xi,X2) in ()2.19|1 . we obtain 



g[xi,X2} 



1 Z"^'' 

— / {r^^ sin 9 - r^^ cos 9) d9, 
47r Jo 



(3.10) 



where r and p are given by ()2.1()j) and 



r(r, p, ^) = /(r, p, /^"^ -h'^ + 2r + /^'"'^ -/i^ - 2/ 



TT 



TT 



(3.11) 



For the numerical calculation of the Hilbert transform we write 



hip, 9) 



' Kp, 0) , , , f hp', 9) - fip, 9) ^ , 



1 P - P 



-dp' + 



P' - P 



dp' 



fip,0)\n 



1-P 
1 + P 



n-l 



f"^^ S.ip',9)-f{p,9) , 
-dp. 



If p = Pi or p = pi+i the integral in the right-hand side of ()3.12|1 can be written 

S,{p',9)-S,{p,f^^ 



(3.12) 



p'-p 



-dp' 



Thus, after some calculations, we obtain 



cip — —Ji + Ji+l 

p - p 



^36 ^^^'^ ~ ^PiPi+^ - 5Pi+i^ - 3(pi - 5pi+i)p - 6p^) 
+^ (5pi^ + 5piPi+i - 4pi+i^ - 3(5pi - pi+i)p + 6p^) f-'^^. 
If p ^ Pi and p 7^ p.j+i the integral in the right-hand side of ()3.12|) can be written 



(3.13) 



p'-p 

and after some calculation we obtain 

n-l 



-dp'-/(p,e) In 



Pi+i - P 



Pi- P 



h{p,9) = 



i=l 
1 



1 



Pi - Pi+l 



In 



Pi+i - P 



Pi- P 



{pi+i - P)fi- {pi - P)fi+i 



6 



{pi - P){pi+i - P) {pt - 2pi+i + p)f'i + (2pi - Pi+i - p)/, 



i+l 



(3.14) 



where Fj is the right-hand side of ()3.13p . 
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In order to calculate numerically J(r, p, 0) for any Xi, X2, 0, we use relations p.ip and 
jnnb). Thus 



f{xi,X2) = -^^ I hp{x2 cost - xi sint, t)dt, 







and consequently 

1 Z"^'' 

F{T,p,e) = hp{Tsm{e-t)+ pcos{e-t),t)dt, (3.15) 

47r^ Jo 

where r and p are given from ()2.1()j) and hp from ()3.4j) . We can now calculate F{t,p,6) 
following the procedure outlined in the previous section. We then calculate J(r, p, 0) using 
relation ()3.5p if r > 0, alternatively the relation 



I{t.P.O) = exp 



f{p,0)- r F(r',p,^)dr' 



(3.16) 



if r < 0. For the numerical calculation of the integrals appearing in ()3.5|1 and ()3.16|1 we 
use the Gauss-Legendre quadrature with two functional evaluations at every step, i.e. 

/ F(r', p, e)dT' ^ w,F{n, p, e) + W2F{t2, p, 9), 

J a 

where the abscissas ri, T2 and the weights wi, W2 are given by 

Ti = a + (/^ - a) I - - — I , Ta = a + (/:/ - a) I - + — I , wx = W2 = -(/? - a). 

We also notice that we have tried subdivision of the interval (a,/3) into several intervals 
and the improvement is very minor. Therefore we use just one interval, i.e. two function 
evaluations per quadrature, since the major increase in running time of the program 
implicit in using panel quadrature is not justified by the modest improvement in accuracy. 

For the numerical calculation of the integrals in ()3.10p and ()3.15p we use again for- 
mula ()3.2|) . resulting in spectral convergence. For the numerical calculation of the partial 
derivatives r^^ and r^^ in ()3.1()j) we use the forward difference scheme 

/'(a;) ~ -3/(a:)+4/ (x + Ax) - f{x + 2Ax) 



2Aa; 

for the first half of the interval [—1,1], and the backward difference scheme 

_ 3/(x)-4/(x-Ax) + /(a;-2Ax) 
^ ^""^ ^ 2Ax 

for the second half. 

Thus, for the numerical calculation of g{xi,X2) from the data f{p,0) and gf{p,6) 
we apply the following procedure: First we calculate the second derivatives //', using 
subroutine spline. Consequently, we calculate h{p,9) using ()3.12p and ()3.13p for all 
given p and 9. We note that if \pi\ = 1, then, since we have assumed compact support, 
f{p,9) = 0, thus the first term in ()3.12j) is absent. We then calculate f'^^^{p,9) and 
f^^'^{p,9) using ()3.6p . as well as f^{p,9) and f'^{p,9) using ()3.8|) (at this stage we use the 
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second data function gf). Finally we calculate, again using spline, the second derivatives 
for the natural cubic spline interpolation of the functions /^(p, 0) and f^{p,6). 

Having calculated all the necessary second derivatives we now proceed as follows: First 
we calculate f{p,9) for any xi, X2 (and 9) using ()2.1()j) and (jH.Hj) . For this purpose we 
have used subroutine splint from Numerical Recipes. Consequently we calculate h{p,9) 
using dnm. Then we calculate /^'""(p, 9) and /'"""{p, 9) using dHUI), /'(p, 9) and /^(p, 9) 
using splint and finally /i^(p, 9) and /i**(p, 9) using relations similar to ()3.14|1 . These last 
six functions are used in (jH.llll . We then calculate I{r,p,9) as described earlier. Finally 
we calculate r{T,p,9) using (jH.llj) and consequently g{xi,X2) using (jH.lOj) . 

4 Numerical Tests 

The 9 points are equally spaced in [0,27r], while the p points are equally spaced in 
[—1, 1]. The density plots presented below were drawn by using Mathematica The 
dark color represents zero (or negative) values while the white color represents the maxi- 
mum value of the original (or reconstructed) function. 

First we tested the PET algorithm for the three different phantoms shown in Figures 
El Figures (a) and (b) were taken from [2H] and [SHj , respectively. These figures depict the 
attenuation coefficient for a function f{xi, X2) modelling a section of a human thorax. The 
small circles represent bones and the larger ellipses the lungs. Figure (c) is the well known 
Shepp-Logan phantom, which provides a model of a head section. All these phantoms 
consist of different ellipses with various densities. 

Using the Radon transform ()1.4|1 . we computed the data function /(p, 9) for 200 points 
for 9 and 100 points for p. This computation was carried out by using Mathematica. We 
then used these data in the numerical algorithm to reevaluate f{xi,X2)- Furthermore, in 
order to remove the effect of the Gibbs-Wilbraham phenomenon, we applied an averaging 
filter as follows: We first found the maximum value (max) of f{xi, X2) in the reconstructed 
image. We then set to zero those values of /(xi, X2) which were less than ^ max. Finally 
we applied the averaging filter with averaging parameter a = 0.005. This filtering proce- 
dure was applied five times, with the additional ehmination of those values of f{xi,X2) 
which were less than ^ max at the end of the procedure. In Figures |3 and El we present 
the results before and after the filtering procedure, respectively. The reconstruction took 
place in a 500 x 500 grid. 

We then tested the SPECT algorithm for the three different phantoms shown in Figures 
El Figures (a) and (b) were taken from |2H|- In these cases the function f{xi,X2) is given 
by Figure El^a). Figure (c) was taken from [20]. The white ring represents the distribution 
of the radiopharmaceutical at the myocardium. In this case the function /(xi, X2) is given 
by Figure E^b). 

By using the Radon transform p.4j) . and the attenuated Radon transform ()1.5|1 . we 
computed the data functions f{p,9) and gf{p,9) for 200 values of 9 and 100 points of 
p (again using Mathematica). We consequently used these data in our program to re- 
evaluate g{xi,X2)- In order to remove the effect of the Gibbs-Wilbraham phenomenon, 
a median filter was used, with the additional elimination of those values of g{xi,X2) 
which were less than ^ max before and after the application of the filter. The results 
are shown in Figures 13 and |H1 before and after the filtering procedure respectively. The 
reconstruction took place in a 140 x 140 grid. 
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(b) 

Figure 3: Test phantoms for the PET algorithm. 




(a) (b) (c) 

Figure 5: The reconstruction of the phantoms of Figures El after the filtering procedure. 



For the above phantoms it seems that even a rough estimation of -F(t, p, 6) is suffi- 
cient for an accurate reconstruction. This means that, in order to compute numerically 
F{t,p,9) using ()3.15|) . it is sufficient to use ten equally spaced points for t, rather than 
200. This reduces considerably the reconstruction time. 
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Figure 6: Test phantoms for the SPECT algorithm. In Figures (a) and (b) the function 
f{xi,X2) is given by Figure Efa), while in Figure (c) the function f{xi,X2) is given by 
Figure ini^b). 




laj (b) [c) 

Figure 7: The reconstruction of the phantoms of Figures El before the filtering procedure. 




(a) (b) [c) 

Figure 8: The reconstruction of the phantoms of Figures IHl after the filtering procedure. 

for useful suggestions. 

References 

[1] D. Koh, G.J.R. Cook, J.E. Husband, New Horizons in Oncologic Imaging (editorial), 
N. Engl. J. Med. 348, 2487 (2003). 

[2] J. Jonides et al.. Verbal and Spatial Working Memory in Humans, Psychol. Learn. 
Motiv. 35, 43 (1996). 



15 



[3] G.S. Mark et al., Understanding Emotional Prosody Activates Right Hemisphere 
Regions, Arch. NeuroL 53, 665 (1996). 

[4] S. Vorstrup, O.B. Paulson, N.A. Lassen, Cerebral Blood Flow in Acute and Chronic 
Ischemic Stroke using Xenon-133 Inhalation Tomography, Acta Neurol. Scand. 74, 
439 (1986). 

[5] M. Lauritzen, J. Olesen, Regional Cerebral Blood Flow During Migraine Attacks by 
Xenon-133 Inhalation and Emission Tomography, Brain 107, 447 (1984). 

[6] B.I. Lee et al., HIPDM-SPECT in Patients with Medically Intractable Complex 
Partial Seizures: Ictal study. Arch. Neurol. 45, 397 (1988). 

[7] J.L. Tyler, T.N. Byme, Neoplastic Disorders, in Clinical Brain Imaging: Principles 
and Applications, eds. J.C. Mazziotta, S. Oilman, p 166, Philadelphia: F.A. Davis 
(1992). 

[8] J.C. Mazziotta, Movement Disorders, in Clinical Brain Imaging: Principles and Ap- 
plications, eds. J.C. Mazziotta, S. Oilman, p 244, Philadelphia: F.A. Davis (1992). 

[9] S. Minoshima et al., A Diagnostic Approach in Alzheimer's Disease Using Three- 
Dimensional Stereotactic Surface Projections of Fluorine-18-FDO PET, J. Nucl. 
Med. 36, 1238 (1995). 

[10] L. Junck et al., PET Imaging of Human Oliomas with Ligands for the Peripheral 
Benzodiazepine Binding Site, Ann. Neurol. 26, 752 (1989). 

[11] J.C. Mazziotta et al.. Reduced Cerebral Olucose Metabohsm in Asymptomatic Sub- 
jects at Risk for Huntington's Disease, N. Engl. J. Med. 316, 357 (1987). 

[12] N.C. Andreasen, Linking Mind and Brain in the Study of Mental Illnesses: A Project 
for a Scientific Psychopathology, Sci. 275, 1586 (1997). 

[13] E.M. Reiman et al., Neuroanatomical Correlates of Anticipatory Anxiety, Sci. 243, 
1071 (1989). 

[14] J.O. Tjuvajev et al., A Oeneral Approach to the Non-Invasive Imaging of Transgenes 
using Cis-Linked Herpes Simplex Virus Thymidine Kinase, Neoplasia 1, 315 (1999). 

[15] Y. Yu et al.. Quantification of Target Oene Expression by Imaging Reporter Oene 
Expression in Living Animals, Nature Med. 6, 933 (2000). 

[16] L.A. Green et al.. Indirect Monitoring of Endogenous Gene Expression by Positron 
Emission Tomography (PET) Imaging of Reporter Gene Expression in Transgenic 
Mice, Mol. Imaging Biol. 4, 71 (2002). 

[17] M. Doubrovin et al.. Imaging Transcriptional Regulation of p53-Dependent Genes 
with Positron Emission Tomography in Vivo, Proc. Natl Acad. Sci. USA, 98, 9300 
(2001). 

[18] D. Lardinois et al.. Staging of Non-Small-Cell Lung Cancer with Integrated 
Positron-Emission Tomography and Computed Tomography, N. Engl. J. Med. 348, 
2500 (2003). 



16 



[19] D. Ost, A.M. Fein, S.H. Feinsilver, The Solitary Pulmonary Nodule, N. Engl. J. Med. 
348, 2535 (2003). 

[20] B.F. Hutton, Cardiac Single-Photon Emission Tomography: Is Attenuation Correc- 
tion Enough? (invited editorial), Eur. J. Nucl. Med. 24, 713 (1997). 

[21] F.J.T. Wackers, Attenuation Correction, or the Emperor's new Clothes? (editorial), 
J. Nucl. Med. 40, 1310 (1999). 

[22] F.M. Bengel et al.. Effect of Sympathetic Reinnervation on Cardiac Performance 
after Heart Transplantation, N. Engl. J. Med. 345, 731 (2001). 

[23] F. Natterer, The Mathematics of Computerized Tomography, Wiley, New York (1986). 

[24] A.S. Fokas, R.G. Novikov, Discrete Analogues of 9-Equations and of Radon Trans- 
form, C. R. Acad. Sci. Paris Ser. I. Math. 313, 75 (1991). 

[25] R.G. Novikov, An Inversion Formula for the Attenuated X-ray Transformation, Ark. 
Mat. 40, 145 (2002). 

[26] L.A. Shepp, B.F. Logan, The Fourier Reconstruction of a Head Section, IEEE Trans. 
Nucl. Sci. 21, 21 (1974). 

[27] F. Natterer, Inversion of the Attenuated Radon Transform, Inv. Prob. 17, 113 (2001). 

[28] L.A. Kunyansky, A New SPECT Reconstruction Algorithm Based on the Novikov 
Exphcit Inversion Formula, Inv. Prob. 17, 293 (2001). 

[29] J. P. Guillement, F. Jaubcrtcau, L. Kunyansky, R. Novikov, R. Trebossen, On Single 
Photon Emission Computed Tomography Imaging based on an Exact Formula for 
the Nonuniform Attenuation Correction, Inv. Prob. 18, Lll (2002). 

[30] J. P. Guillement, R.G. Novikov, A Noise Property Analysis of Single-Photon Emission 
Computed Tomography Data, Inv. Prob. 20, 175 (2004). 

[31] T. Hcbcrt, R. Leahy, M. Singh, Fast MLE for SPECT using an Intermediate Polar 
Representation and a Stopping Criterion, IEEE Trans. Nucl. Sci. 35, 615 (1988). 

[32] Z. Liang, H. Hart, Bayesian Reconstruction in Emission Computed Tomography, 
IEEE Trans. Nucl. Sci. 35, 788 (1988). 

[33] J. Nuyts, J. A. Fessler, A Penalized-Likelihood Image Reconstruction Method for 
Emission Tomography, compared to Post-Smoothed Maximum-Likelihood with 
Mached Spatial Resolution, IEEE Trans. Med. Imag. 22, 1042 (2003). 

[34] A.S. Fokas, I.M. Gel'fand, Intcgrability of Linear and Nonlinear Evolution Equations, 
and the Associated Nonlinear Fourier Transforms, Lett. Math. Phys. 32, 189 (1994). 

[35] M.J. Ablowitz, A.S. Fokas, Introduction and Applications of Complex Variables, Cam- 
bridge University Press (1997). 

[36] B. Fornberg, A Practical Guide to Pseudospectral Methods, Cambridge University 
Press (1996). 



17 



[37] W.H. Press, S.A. Teukolsky, W.T. Vcttcrling, B.P. Flannery, Numerical Recipes in 
Fortran. The Art of Scientific Computing (2nd edition), Cambridge University Press 
(1992). 

[38] S. Wolfram, The Mathematica Book (4th edition), Cambridge University Press (1999). 



18 



